# load packages
library(tidyverse)
library(flextable)
library(plm)
library(sjPlot)

# import the data
data <- readRDS("data.RDS")

# Summary Statistics
summary2 <- function(x) {
  N <- map_int(x, function(x) sum(!is.na(x)))
  Mean <- map_dbl(x, mean, na.rm = TRUE)
  S.D. <- map_dbl(x, sd, na.rm = TRUE)
  Min. <- map_dbl(x, min, na.rm = TRUE)
  Max. <- map_dbl(x, max, na.rm = TRUE)
  tibble(Variable = colnames(x), N, Mean, S.D., Min., Max.)
}

sum_stats <- data %>%
  select(fariss.mean, ai_shaming, us_shaming, ln_pop, ln_gdppc, mediascore, democracy, civwar) %>%
  summary2() %>% flextable() %>%
  colformat_double(j = 3:6, digits = 2)

# Table C1
sum_stats

# regression analysis
model1 <- lm(fariss_next ~ ai_shaming + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

model2 <- lm(fariss_next ~ us_shaming + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

model3 <- lm(fariss_next ~ ai_shaming + us_shaming + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

pdata <- pdata.frame(data, index = c("country_cow_code", "year"))

model4 <- plm(fariss_next ~ ai_shaming + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

model5 <- plm(fariss_next ~ us_shaming + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

model6 <- plm(fariss_next ~ ai_shaming + us_shaming + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

# Figure 4
plot_models(model1, model2, model3, rm.terms = c("fariss.mean", "ln_gdppc", "ln_pop", "mediascore", "democracy", "civwar"), m.labels = c(paste("Model", 1:3))) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylim(-.04, .04) +
  theme_bw() +
  theme(legend.position = "top", legend.title = element_blank())
ggsave("Figures/fig4.jpg", width = 5, height = 4)

# Table E1
tab_model(model1, model2, model3, model4, model5, model6, collapse.se = TRUE, show.ci = FALSE, p.style = "stars",
          dv.labels = str_c("Model ", 1:6), digits = 3)

# Robustness Check
## Alternative Measurement
model1 <- lm(fariss_next ~ ai_shaming2 + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

model2 <- lm(fariss_next ~ us_shaming2 + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

model3 <- lm(fariss_next ~ ai_shaming2 + us_shaming2 + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

pdata <- pdata.frame(data, index = c("country_cow_code", "year"))

model4 <- plm(fariss_next ~ ai_shaming2 + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

model5 <- plm(fariss_next ~ us_shaming2 + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

model6 <- plm(fariss_next ~ ai_shaming2 + us_shaming2 + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

### Table E2
tab_model(model1, model2, model3, model4, model5, model6, collapse.se = TRUE, show.ci = FALSE, p.style = "stars",
          dv.labels = str_c("Model ", 1:6), digits = 3)

## Alternative DV
### generate the average score of future 10 years
fariss <- read_csv("HumanRightsProtectionScores_v4.01.csv")

data$fariss_average <- NA

for (i in 1:nrow(data)) {
  cat(data$country_cow_code[i], data$year[i], "\n")
  data$fariss_average[i] <- mean(fariss$theta_mean[fariss$YEAR > data$year[i] & fariss$YEAR < (data$year[i] + 11) & fariss$COW == data$country_cow_code[i]])
}

model1 <- lm(fariss_average ~ ai_shaming + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

model2 <- lm(fariss_average ~ us_shaming + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

model3 <- lm(fariss_average ~ ai_shaming + us_shaming + fariss.mean + ln_pop + ln_gdppc + mediascore + democracy + civwar, data = data)

pdata <- pdata.frame(data, index = c("country_cow_code", "year"))

model4 <- plm(fariss_average ~ ai_shaming + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

model5 <- plm(fariss_average ~ us_shaming + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

model6 <- plm(fariss_average ~ ai_shaming + us_shaming + fariss.mean + ln_gdppc + ln_pop + mediascore + civwar + democracy, model = "within", 
              effect = "individual", data = pdata)

## Table E3
tab_model(model1, model2, model3, model4, model5, model6, collapse.se = TRUE, show.ci = FALSE, p.style = "stars",
          dv.labels = str_c("Model ", 1:6), digits = 3)